Particle currents and the distribution of terrace sizes in unstable epitaxial growth 
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A solid-on-solid model of epitaxial growth in 1 + 1 dimensions is investigated in which slope depen- 
dent upward and downward particle currents compete on the surface. The microscopic mechanisms 
which give rise to these currents are the smoothening incorporation of particles upon deposition and 
an Ehrlich-Schwoebel barrier which hinders inter-layer transport at step edges. We calculate the 
distribution of terrace sizes and the resulting currents on a stepped surface with a given inclina- 
tion angle. The cancellation of the competing effects leads to the selection of a stable magic slope. 
Simulation results are in very good agreement with the theoretical findings. 
PACS numbers: 81.10.Aj, 05.70.Ln, 68.55.-a 



Epitaxial growth has become a standard method for 
the production of high-quality crystals and films, needed 
for e.g. semiconductor devices. An overview of experi- 
mental techniques can be found in |5J], for instance. Sig- 
nificant effort has been devoted to a theoretical under- 
standing of the many morphologies and scaling behaviors 
that can be observed in epitaxial growth, see e.g. || for 
a review of theoretical approaches. Here we address the 
frequently observed phenomenon of mounds in unstable 
growth which has attracted considerable interest, see e.g. 
[^-^|. Specifically, we consider situations in which com- 
peting smoothening and steepening effects control the 
surface morphology and lead to the selection of a stable 
slope in the system. 

We discuss potential microscopic mechanisms which re- 
sult in the emergence of mounds and slope selection in 
the frame of a discrete (1 + l)-dimensional model. The 
net particle currents as well as the distribution of terrace 
sizes on a surface of a given slope can be worked out and 
this allows then to evaluate the magic slope as well as the 
complete statistical properties of the emerging surface. 
The analysis complements previous theoretical investiga- 
tions which address the mean terrace size only or neglect 



fluctuations explicitly 
Frank (BCF) theory ~ 



in the spirit of Burton Cabrera 
LfJ . We demonstrate that the full 
distribution of terrace sizes carries relevant information 
that should be taken into account. Our results suggest, 
for example, that it should be possible to identify relevant 
microscopic mechanisms from experimental data. 

Our (1 + l)-dimensional model obeys the solid-on- 
solid (SOS) restriction, i.e. the surface can be described 
by an integer array of height variables h k . Single 
particles are deposited at randomly chosen sites k £ 
{1, 2, . . . , L}. Upon arrival, an incorporation process 
moves a particle to the lowest available site within a 
neighborhood of ±i? lattice constants, i.e. the site j 



.{h 



k-R, 



, h k , . . . h k+R ). In case of a 



with hj = 

tie, the site closest to the deposition is chosen. Only 



if this is still ambiguous, an additional random selec- 
tion is performed. Such a smoothening mechanism, in 
absence of further effects, is commonly associated with 
the Edwards-Wilkinson universality class of growth [||. 
The parameter R — 0(1) (in lattice constants) is termed 
the incorporation radius and sets the typical length scale 
of the process. Various interpretations of incorporation 
have been considered, including downward funneling on 
non-trivial lattices and knock-out -processes due to the 
momentum of incoming particles, see e.g. pj] , p^ |. 

A particle which is, after deposition and possibly in- 
corporation, not yet bound to a lateral neighbor diffuses 
on the surface by performing a random walk (RW) until 
it reaches an additional binding partner and becomes im- 
mobile or until it collides with another moving adatom 
and forms an island nucleus. In a density of diffusing 
particles, this nucleation process would result in a typical 
collision free path Id- In the case of irreversible aggre- 
gation on a flat substrate, and if islands of two or more 
adatoms are considered immobile, it has been shown that 
Id en (d/f) 1 ^ 4 in (1 + 1) dimensions j^|l]|. Here, d is 
the diffusion constant and / the incoming flux, with all 
lengths in dimensionless lattice constants. 

In J|-§] nucleation is represented in an effective sin- 
gle particle picture and Id fixes the typical distance of 
island nuclei in the first layers on a flat substrate. After 
the formation of mounds, terrace sizes are much smaller 
than Id, typically. Hence we will completely neglect nu- 
cleation on the stepped surface which is justified for small 
incoming flux or fast diffusion, respectively. 

Here, the RW ends whenever the particle sticks ir- 
reversibly to a lateral neighbor, i.e. when it reaches 
a terrace step. Attachment can occur from below and 
above, in principle, but this symmetry is broken due to 
the so-called Ehrlich-Schwoebel (ES) effect [§: An ad- 
ditional energy barrier E es at step edges hinders down- 
ward moves of diffusing adatoms as this would involve 
loosely bound intermediate positions. Our model takes 
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the ES effect into account by assigning a probability 
p es oc exp [— E es /(kT)] to downward moves. 

The ES effect results in an uphill current of adatoms 
because particles will stick to upper terraces preferably 
for any p es < 1. On the other hand, incorporation con- 
stitutes a downhill current. Both effects are slope depen- 
dent and their cancellation gives rise to the formation of 
mounds with a well-defined inclination angle. Once the 
structures have built up the magic slope, a coarsening 
process begins which decreases the number of mounds, 
see H H for details. 

We will first work out the distribution of terrace sizes 
which emerges in our model on surfaces with a given, 
fixed inclination. Further, we will calculate the mean 
displacement of a particle's final position from its depo- 
sition site. The latter corresponds to the net particle 
current on the growing surface and its zero as a function 
of the inclination angle determines the stable slope. 

We consider a triple of terraces on an inclined sur- 
face with a central terrace c of width l c (lattice sites 
j = 1,2,... l c ) and its neighbors t to the left (r to the 
right) with size le (l r ). Without loss of generality, we as- 
sume that the surface height decreases to the right. For 
a particle deposited on a site k of terrace c with l c > R 
we have to distinguish the following cases: 

(a) k > l c — R: the incorporation process places the 
particle at its final position at site l c + 1, attaching 
to the lower terrace end. 

(b) 1 < fc < l c — R: the particle performs a random 
walk until it reaches one of the trap sites j = 1 or 
j = l c + l where it comes to rest. This includes de- 
position at site k — 1 without subsequent diffusion. 

A diffusing particle located at a site j with 2 < j < l c ~ 1 
moves to one of the neighboring positions j ± 1 with equal 
probability 1/2. The asymmetry of the RW (b) is due to 
the ES-barrier present for jumps from site j — l c . We 
denote with Pes/2 the probability for a downward move 
to site l c + 1. With probability (1 — Pes) 1 2 the move is 
rejected and the particle remains at j = l c for the next 
time step, whereas with probability 1/2 it jumps to l c — 1. 

A straightforward exercise yields the probability for a 
RW initiated at site k to end in j = l c + l by downward 
diffusion: q(k,l c ) = ((k-l)p es )/ ( l + (l c -l)p es ) which 
obviously satisfies q(l,l c ) = for all p es . For similar 
problems of this type see for instance Q . 

Hence, the total probability d(l c ) for a deposition event 
to occur on terrace c with subsequent downward diffusion 
is d(l c ) = A(l c )/L with 

y ' fr[ q{ ' ; 2 + 2(l c -i) Pes 1 ' 

if l c > R and A(l c ) = else. The second case accounts 
for the fact that any particle directly deposited onto a 
terrace of width l c < R will be incorporated without 
performing diffusion. 



The quantity A(7 C ) can be interpreted as the effective 
number of deposition sites which contribute to downward 
diffusion from terrace c. The prefactor 1/L of al(l c ) is 
simply the constant probability for deposition on any of 
the sites in the system. 

Now we can work out the probability w(l c ^l c — 1) for 
the central terrace c to be shortened by the next deposi- 
tion event. For l c > R one finds 

w(l c ^l c -l) = (R+[l c -R-A(l c )] + A{k))/L 

= (l c -A(l c )+A(k))/L. (2) 

The first contribution, R, represents deposition on any 
of the R sites left of terrace c. Note that the outcome 
of the subsequent incorporation process is completely in- 
dependent of the surface configuration, in particular of 
the left neighbor terrace width It. The second term [ . . . ] 
accounts for deposition on c with final attachment to the 
upper terrace. Finally A(lg)/L is the probability for the 
shortening of c through diffusion from terrace t. 

If < l c < R, only two processes can shorten the 
central terrace: incorporation from exactly l c sites left of 
c and downward diffusion from terrace £. One obtains 
w(l c —>lc— 1) = {lc+A(le))/L and we finally observe that 
Eq. (||) is valid for any l c > 1 since A(Z C ) = for all 
l c < R. Obviously, w(l c — > l c — 1) = if l c = already. 

Since l c can only increase at the cost of shortening l r 
at the same time, one obtains immediately the result 

^ lc+1) = ^ L (l r -A(l,) + A(l c )) ifj : >0 (3) 

We proceed by assuming that in a population of terraces 
the distribution of their sizes factorizes, i.e. that a single, 
identical p(l) is sufficient to describe their statistics. For 
the limiting case of an infinite ES barrier (p es — 0) the 
evolution of a terrace is independent of the entire config- 
uration left of it, and the above property can be shown 
to hold true. In general, neighbor terraces clearly inter- 
act. Nevertheless, our simulations show that the assump- 
tion of identically distributed, independent terrace sizes 
yields excellent approximations, at the very least. Figure 
|l| shows, for instance, that the correlation coefficient of 
neighboring terrace sizes in a system with L — 1000 and 
p es — 0.2 vanishes within error bars. 

The analysis simplifies significantly if terrace sizes are 
considered to be uncorrelated. The above expressions 
(||,|]) were obtained for a given triple of terraces { £, c, r }. 
By averaging w(l c — > / c ±l) over p{h) and p(l r ), respec- 
tively, one obtains the mean probabilities w(l c ~^l c ±l): 

, i(l c -A(l c ) + (A)) for l c >0 

w(i c ^L c for ? c = [ > 

w(l c ^l c + l)= (<J>-<A) + [l-p(0)]A(Q)/L (5) 

where the r.h.s. involve only the width of the consid- 
ered terrace itself, the frequency p(0) of vanishing ter- 
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race sizes and the mean values (l) 
(A>=E,°lo AO>(j),see Eq. f). 

The evolution of terraces according to (|4]||) produces 
a stationary distribution p(l), if p(l + l)w(l + l — > /) = 
p(Z)u>(Z —►£+!), hence 



p(i + i) = p(0 



(0-(A) + [l-p(0)]A(Q 
1 + 1- A(Z + 1) + (A) 



(6) 



for I > 0. This relation is implicit, since all have to 
be known for the evaluation of the averages on the r.h.s. 
In the particular case of an infinite ES barrier, Eq. (|^) 
reads p(l + l) = p(l) (I) /{I + 1) which is satisfied by the 
Poissonian p(l) = X 1 e~ x /l\ with mean (I) = A. 

In order to obtain the stationary p(l) for p es > on 
a surface with a given mean terrace size A, we replace 
(I) with A and (A) with an adjustable parameter D in 
Eq. (Q). The quantities p(0) and D are determined such 
that ^2 t p(l) — 1 is satisfied and (/) = A is reproduced 
self-consistently. Note that then, by construction, (JF^) 
guarantees (A) — D as well. In the numerical treatment, 
sums are truncated at a value l m ax, with the resulting 
p(Z max ) small enough to justify the truncation a posteri- 
ori. 

A particle deposited at, say, lattice site i will become 
immobile at a final position j ^ i after incorporation 
and diffusion, in general. The expected displacement 5 = 
{j — i) depends on the distribution of terrace sizes in the 
system. Taking into account all possible displacement 
processes and their corresponding probabilities one finds 



5=\r{R 



i)-i|>(0[(fl-0(fl-i + i)]+ (7) 

1=0 



1. 
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£ Pit) -l(l-R-A(l)) + -l(l + l) - -R(R+l) 

l=R+l ^ 

where the p(l) obtained from (Q) have to be inserted for a 
given (I) — A. Here, the first line corresponds to the ex- 
pected (positive, downward) effect of incorporation and 
the second represents the total (negative, upward) con- 
tribution of diffusion. In the limiting case p es = Eq. 
(0) reduces to 6 = R(l) — (I) 2 /2, exploiting the fact that 
A(Z) = in this case and (J 2 ) — (l) 2 = (I) for the Poisson 
distribution. 

Figure [j] shows the result of Monte Carlo simulations 
of the growth process in a system of size L = 1000 for 
the model with R — 2 and different values of p es , where 
boundary conditions were used to fix (/) = A. We have 
displayed the particle current 5/X on the surface as a 
function of A. Note that S in Eq. (0) was obtained for the 
normalization YliPQ) = 1 which corresponds to a fixed 
number of terraces. In systems with a fixed number L of 
lattice sites, an additional factor has to be introduced as 
the number of terraces grows like 1/A. 

For very steep surfaces, A — > 0, the displacement ap- 
proaches the limiting value R, representing the fact that 




FIG. 1. The mean particle displacement S/X vs. A = (I), 
cf. Eq. (H), shown for the model with R — 2 and different 
ES barriers. Symbols represent the result of simulations with 
L = 1000 for Pes = 0(«), 0.1(B), 0.2(T), 0.4(A). Averages 
were performed over 100 runs and standard error bars would 
be smaller than the symbols, where not shown. In addition, 
the correlation coefficient p = ({lalb} — (l) 2 ) / ((l 2 ) — (l) 2 ) for 
neighboring terraces a and b is displayed for the case p es = 0.2; 
crosses and error bars correspond to 10 • p. 



every deposited particle is shifted by R lattice constant 
in the incorporation and then comes to rest. In the limit 
of vanishing slope, our model yields a diverging negative 
"upward" current. This is an artifact of completely ne- 
glecting nucleation, which inevitably becomes important 
as A — > oo and imposes a maximal displacement on the 
order of In- 

On mounded surfaces, bottom and top terraces limit 
the extension of inclined flanks. Any slope that results 
in a net uphill current according to Eq. (|?]) will steepen 
this portion of the surface and vice versa. Accordingly, 
a mean terrace width A Q will be stabilized which corre- 
sponds to the zero of (5(A). In the presence of an infinite 
ES-barrier we find the exact relation A = 2R. A naive 
and not quite correct argument was used in || to ob- 
tain the same result. It is instructive to check that the 
magic slope cannot be obtained from the condition that 
the mean displacement vanishes on a particular terrace of 
size A. This would correspond to setting p{l) — 5 l ^ in Eq. 
(0) and gives results analogous to the BCF-like treat- 
ment in H which does not account for fluctuating terrace 
sizes. For p es = one obtains, e.g., A = 2i?+l = A + l. 

Fig. H shows the frequency of terrace sizes as observed 
in two different systems which both stabilize the mean 
A = 6. In one we have set p es = and R = 3, the 
second example corresponds to R = 2 and p es = 0.258. 
Computer simulations show excellent agreement. Sys- 
tems with a very pronounced ES-effect produce a narrow 
distribution with a very low frequency p(0) of step bunch- 
ing, i.e. zero terrace sizes. As a limiting case one finds 
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for infinite ES-barrier (j> es = 0). 
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FIG. 2. Frequency p(l) of terrace sizes in two cases with 
(?) = A = 6. The dashed line shows the theoretical prediction 
for R — 3, Pes = 0, triangles represent simulations (as in Fig. 
[l]); the solid line and bullets correspond to R = 2,p ea = 0.258. 
The dot-dashed curve displays a geometric distribution with 
{1} = 6 for comparison. 
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FIG. 3. The selected mean A (solid lines) and 
a 2 = (l 2 ) — \ 2 (dashed) vs. p es . Pairs of curves correspond 
to (from below) R — 1,2, 3. Note that for p es = we find 
A = a = 2R. The symbols represent two choices of (p es , R) 
which result in A = 6, cf. Fig. |j| 



On the contrary, step bunching is observed with a much 
larger frequency in cases with a weaker ES-effect where 
the distribution p(l) is much broader. Fig. ^ displays 
A„ and the variance a 2 of the terrace size distribution 
as functions of p es . Note that a 2 grows drastically with 
increasing p es , indicating significant deviations from the 
Poissonian for infinite ES-effect. 

The analysis of experimental data is frequently based 
on the simple assumption of random, non-interacting ter- 
race sizes on vicinal surfaces. This leads to the geometric 
distribution p g (l) = (1 — 1/A)' _1 /A with (I) = A, see e.g. 
|15| for a discussion. Note that p g (l) differs significantly 
from the type of statistics that we find in our model, cf. 
Fig. |^. In particular, step bunching is much more fre- 
quent in this simple picture: p g (0) — (A— 1) . 

In summary we have presented a microscopic model of 
unstable epitaxial growth in which it is possible to derive 
the net particle currents on surfaces of a given inclination. 
For the first time it is possible to work out the full distri- 
bution of terrace sizes in such a system. Further, we were 
able to calculate the stable mean terrace size and the cor- 
responding statistical properties of the surface. We have 
restricted ourselves to the analysis of (1 + l)-dimensional 
growth in this work. However, our results should carry 
over to a more realistic (2 + 1) -dimensional picture to 
a large extent, whenever terrace edges do not meander 
significantly. 

Our findings allow for a qualitative interpretation of 
experimental results in systems which display slope se- 
lection: frequent step bunching and a broad distribution 
hint at a relatively weak ES-barrier. Narrow distribu- 
tions with little or no step bunching indicate that a sig- 
nificant ES-effect is present but is compensated for by 
smoothening effects like downhill tunneling. 

Extensions of this work will concern desorption and its 



influence on the growth process. Preliminary results indi- 
cate that a significant desorption rate can trigger a tran- 
sition from slope selection to rough growth and we expect 
non-trivial effects in the statistics of terrace sizes. 
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